/*
 * lib.h
 *
 *  Created on: 2010.06.26
 *      Author: Tadas
 */

#ifndef LIB_H_
#define LIB_H_

#include <cmath>

#define EPS 1e-8
#define JMAX 16

using namespace std;

const int lygSkaic = 2;

double stepRK1(double x, double y, double h, double (*func)(double, double)){
	return y + h*func(x, y);
}

double stepRK2(double x, double y, double h, double (*func)(double, double)){
	double k1, k2;
	
	k1 = h * func(x, y);
	k2 = h * func(x + 0.5 * h, y + 0.5 * k1);
	
	return y + k2;
}

double stepRK4(double x, double y, double h, double (*func)(double, double)){
	double k1, k2, k3, k4;
	
	k1 = h * func(x, y);
	k2 = h * func(x + 0.5 * h, y + 0.5 * k1);
	k3 = h * func(x + 0.5 * h, y + 0.5 * k2);
	k4 = h * func(x + h, y + k3);
	
	return y + (k1 + 2.0 * k2 + 2.0 * k3 + k4)/6.0;
}




double xDeriv(double t, double *y){
	return y[1];
}

double yDeriv(double t, double *y){
	return -y[0];
}

double (*funkcija[lygSkaic])(double, double *) = {
		xDeriv,
		yDeriv
};

void stepP(double *y, double t, double h, double (**func)(double, double *)){
	int i;
	double th, tt, t6;
	
	double dym[lygSkaic], dyt[lygSkaic], ytt[lygSkaic], yt[lygSkaic], f[lygSkaic];
	
	tt = 0.5 * h;
	t6 = h/6.0;
	
	th = t + tt;
	
	for (i = 0; i < lygSkaic; i++){
		f[i] = func[i](t, y);
		yt[i] = y[i] + tt*f[i];
	}
	
	for(i = 0; i < lygSkaic; i++){
		dyt[i] = func[i](th, yt);
		ytt[i] = y[i] + tt * dyt[i];
	}
	
	for(i = 0; i < lygSkaic; i++){
		dym[i] = func[i](th, ytt);
		yt[i] = y[i] + h * dym[i];
		dym[i] += dyt[i];
	}
	
	for(i = 0; i < lygSkaic; i++){
		dyt[i] = func[i](t + h, yt);
		y[i] += t6 * (f[i] + 2 * dym[i] + dyt[i]);
	}
}



#endif /* LIB_H_ */
